# TODO: Add comment
# 
# Author: rogb
###############################################################################


setClass("VolatilitySurface",
		contains=c("SituationDate","BasicDateData","Strike")
)

VolatilitySurface <- function(SituationDate,Name,Currency,Date,Strike,Data){
	
	colnames(Data) <- as.character(Strike)
	rownames(Data) <- as.character(as.numeric(Date))
	
	Strike <- sort(Strike)
	Date <- sort(Date)
	
	Data[] <- Data[as.character(as.numeric(Date)),as.character(Strike)]

	new("VolatilitySurface",SituationDate=SituationDate,Name=Name,Currency=Currency,Date=Date,Strike=Strike,Data=Data)
}

# vvv <- matrix(runif(10*6),nrow=10)
# vvv[sample(1:length(vvv),30,replace=F)] <- NA
# v1 <- VolatilitySurface(Sys.Date(),Name="SMI Index",Currency="CHF",Date=Sys.Date()+(1:10)*100,Strike=80+(0:5)*10,vvv)
# v2 <- VolatilitySurface(Sys.Date(),Name="SX5E Index",Currency="CHF",Date=Sys.Date()+(1:10)*100,Strike=80+(0:5)*10,vvv)
# x=1:6, y=1:10 vvv <- t(x %*% t(y))

setMethod("show","VolatilitySurface",function(object){
			tmp <- object@Data
			rownames(tmp) <- getDate.character(object)
			colnames(tmp) <- strike(object)
			cat(name(object)[1]," Volatility Surface\n",situationDate.character(object),"\n",sep="")
			print(tmp)
		})

setMethod("plot",signature=c("VolatilitySurface"),function(x,y,...){
			dimObject <- dim(x)
			TimeToMaturity <- yearFraction(situationDate(x),getDate(x))
			Strike <- strike(x)
			Volatility <- getData(x)
			
			if(dimObject[1]==1){
				main <- paste(name(x)[1]," Volatility Smile\n",as.character(situationDate(x)),sep="")
				plot(Strike,Volatility,lwd=3,col=2,main=main,xlab = "Strike",ylab = "Volatility")
				lines(Strike,Volatility,lwd=2)
			}
			if(dimObject[2]==1){
				main <- paste(name(x)[1]," Volatility Term Structure\n",as.character(situationDate(x)),sep="")
				plot(TimeToMaturity,Volatility,lwd=3,col=2,main=main,xlab = "Time To Maturity",ylab = "Volatility")
				lines(TimeToMaturity,Volatility,lwd=2)
			}
			if(dimObject[1]>1 & dimObject[2]>1){
				main <- paste(name(x)[1]," Volatility Surface\n",as.character(situationDate(x)),sep="")
				persp(TimeToMaturity,Strike,Volatility,theta = -45, phi = 25, expand = 0.5, col = "lightblue", 
						ltheta = 120, shade = 0.75, ticktype = "detailed",
						xlab = "Time To Maturity", ylab = "Strike", zlab = "Volatility",main=main)
			
			}
		})

getVolatilitySurface <- function(Date,Strike,InputVolatilitySurface){
	# radial basis function used for interpolation
	# InputVolatilitySurface <- v2
	# Date <- DateYMD(2011,07,20)
	# Date <- DateYMD(2011,12,14)
	# Date <- as.Date("2011-09-08")
	# Date <- as.Date("2010-11-13")
	# Strike <- 105
	InputSituationDate <- situationDate(InputVolatilitySurface)
	InputName <- name(InputVolatilitySurface)[1]
	InputCurrency <- currency(InputVolatilitySurface)[1]
	InputDate <- getDate(InputVolatilitySurface)
	nInputDate <- length(InputDate)
	
	InputStrike <- strike(InputVolatilitySurface)
	nInputStrike <- length(InputStrike)
	
	Vola <- getData(InputVolatilitySurface)
	
	#dataSet <- data.frame(expand.grid(Date=InputDate,Strike=InputStrike),Vola=as.numeric(Vola))
	Amat <- matrix(NA,ncol=prod(dim(Vola)),nrow=prod(dim(Vola)))
	bVec <- as.numeric(Vola)
	i <- 0
	for(iStrike in InputStrike){
		# iStrike <- InputStrike[1]
		for(iDate in InputDate){
		# iDate <- InputDate[1]
			TimeDifference <- as.numeric(InputDate - iDate)
			TimeDifference <- TimeDifference/(max(TimeDifference)-min(TimeDifference))
			TimeDifference[is.na(TimeDifference)] <- 0
			
			StrikeDifference <- InputStrike - iStrike
			StrikeDifference <- StrikeDifference/(max(StrikeDifference)-min(StrikeDifference))
			StrikeDifference[is.na(StrikeDifference)] <- 0
			
			Dx <- rep(1,nInputDate)%*%t(StrikeDifference)
			Dy <- -TimeDifference%*%t(rep(1,nInputStrike))
			DD <- sqrt(Dx^2 + Dy^2)
			i <- i + 1
			Amat[i,] <- as.numeric(DD)
		}
	}
	Amat <- Amat[!is.na(bVec),!is.na(bVec)]
	bVec <- bVec[!is.na(bVec)]
	
	fit <- solve(Amat,bVec)
	
	dataSet <- expand.grid(Date=as.numeric(Date),Strike=Strike)
	
	result <- apply(dataSet,1,function(x){
		# Date <- as.numeric(Date[1])
		# Strike <- as.numeric(Strike[1])
		Date <- x[1]
		Strike <- x[2]
		TimeDifference <- as.numeric(InputDate - Date)
		TimeDifference <- TimeDifference/(max(TimeDifference)-min(TimeDifference))
		TimeDifference[is.na(TimeDifference)] <- 0
		#TimeDifference <- sapply(Date,function(x,y){tdiff <- as.numeric(y-x);tdiff/(max(tdiff)-min(tdiff))},InputDate)
				
		StrikeDifference <- InputStrike - Strike
		StrikeDifference <- StrikeDifference/(max(StrikeDifference)-min(StrikeDifference))
		StrikeDifference[is.na(StrikeDifference)] <- 0		
		#StrikeDifference <- sapply(Strike,function(x,y){sdiff <- as.numeric(y-x);sdiff/(max(sdiff)-min(sdiff))},InputStrike)
				
		Dx <- rep(1,nInputDate)%*%t(StrikeDifference)
		Dy <- -TimeDifference%*%t(rep(1,nInputStrike))
		DD <-  sqrt(Dx^2 + Dy^2)
		DD[is.na(Vola)] <- NA
		as.numeric(DD[!is.na(DD)]) %*% fit	
	})
	nStrike <- length(Strike)
	VolatilitySurface(InputSituationDate,InputName,InputCurrency,Date,Strike,matrix(result,ncol=nStrike))
}

setMethod("[",signature=c("VolatilitySurface","missing","character"),function(x,i,j,drop){
			Strike <- as.numeric(j)
			Date <- getDate(x)
			getVolatilitySurface(Date,Strike,x)			
		})

setMethod("[",signature=c("VolatilitySurface","Date","missing"),function(x,i,j,drop){			
			Strike <- strike(x)
			Date <- i
			getVolatilitySurface(Date,Strike,x)		
		})

setMethod("[",signature=c("VolatilitySurface","Date","character"),function(x,i,j,drop){			
			Strike <- as.numeric(j)
			Date <- i
			getVolatilitySurface(Date,Strike,x)		
		})

setMethod("[",signature=c("VolatilitySurface","missing","missing"),function(x,i,j,drop){			
			Strike <- strike(x)
			Date <- getDate(x)
			getVolatilitySurface(Date,Strike,x)		
		})

smoothVolatilitySurface <- function(InputVolatilitySurface,...){
	# InputVolatilitySurface <- v1
	
	InputSituationDate <- situationDate(InputVolatilitySurface)
	InputName <- name(InputVolatilitySurface)[1]
	InputCurrency <- currency(InputVolatilitySurface)[1]
	InputDate <- getDate(InputVolatilitySurface)
	nInputDate <- length(InputDate)
	
	InputStrike <- strike(InputVolatilitySurface)
	nInputStrike <- length(InputStrike)
	
	Vola <- getData(InputVolatilitySurface)
	
	TimeDifference <- as.numeric(InputDate - mean(InputDate))
	TimeDifference <- TimeDifference/(max(TimeDifference)-min(TimeDifference))
	TimeDifference[is.na(TimeDifference)] <- 0
			
	StrikeDifference <- InputStrike - mean(InputStrike)
	StrikeDifference <- StrikeDifference/(max(StrikeDifference)-min(StrikeDifference))
	StrikeDifference[is.na(StrikeDifference)] <- 0
	
	dataSet <- data.frame(expand.grid(Date=TimeDifference,Strike=StrikeDifference),Vola=as.numeric(Vola))
	dataSet.clean <- na.omit(dataSet)
	mod <- loess(Vola~Date + Strike,data=dataSet.clean,...)
	vvv <- predict(mod,newdata=dataSet[,-3])
	VolatilitySurface(InputSituationDate,InputName,InputCurrency,InputDate,InputStrike,matrix(vvv,ncol=nInputStrike))
}





